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This article gives an introduction to the multilevel blocking (MLB) approach to both the fermion 
Q and the dynamical sign problem in path-integral Monte Carlo simulations. MLB is able to sub- 

D ■ stantially relieve the sign problem in many situations. Besides an exposition of the method, its 

^0 ' accuracy and several potential pitfalls are discussed, providing guidelines for the proper choice of 

certain MLB parameters. Simulation results are shown for strongly interacting electrons in a 2D 
parabolic quantum dot, the real-time dynamics of several simple model systems, and the dissipative 
two-state dynamics (spin-boson problem). 
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I. INTRODUCTION: THE SIGN PROBLEM 



^ . Quantum Monte Carlo (QMC) techniques are among the most powerful and versatile methods for the computer 
simulation of a large variety of interesting quantum systems encountered in physics, chemistry or biology Q. In 
particular, QMC is capable of delivering numerically exact results. Despite the great potential of this method, there 
are several restrictions and handicaps inherent to all QMC techniques, the perhaps most pressing one being due to 
the sign problem [||. There are various sign problems, namely the fermionic sign problem encountered in equilibrium 
(imaginary-time) simulations of strongly correlated many-fermion systems, and the dynamical sign problem in real- 
time (dynamical) simulations, which already shows up for a single particle. Unfortunately, apart from variational 
i i' or approximate treatments (such as the fixed-node approximation), a completely general and totally satisfactory 
solution to the sign problem in QMC simulations is still lacking. Nevertheless, over the past few years considerable 
and substantial progress has been achieved along several different lines without introducing approximations into the 
, QMC scheme, see, for instance, Refs. [[§-|5). 

■ In these notes focus is put on one specific class of QMC methods called Path-integral Monte Carlo (PIMC). PIMC 
7-H , is based on a discretized path-integral representation of the quantities of interest. The sign problem then arises 
ON ' when different paths that contribute to averages carry different sign (or complex- valued phase). For instance, for 
the fermionic sign problem, as a consequence of exchange, one typically has to deal with determinants, so that 
non-positive-definite fermionic density matrix elements arise. The sign cancellations arising from sampling fermion 
paths then manifest themselves as a signal-to-noise ratio, r\ ~ exp(—7V/3_Eo), that vanishes exponentially with both 
C$ ■ particle number N and inverse temperature (3 — l/fc^T; here Eq is a system-dependent energy scale. The small signal 
surviving the interference of many fermion paths is then inevitably lost in the large background noise of the stochastic 
simulation. Similarly, when studying the real-time dynamics of even a single particle, the quantum-mechanical time 
evolution operator exp(—iHt/K) attaches a complex-valued phase to each quantum path, which in turn gives rise to 
the dynamical sign problem. Again the signal-to-noise ratio will vanish exponentially, r\ ~ exp(— t* /tq), where t* 
, is the maximum real time under study and To a system- and implementation-specific characteristic time scale. The 
1* exponential scaling is typical for the "naive" approach, where one simple uses the absolute value of the complex- valued 
. ! weight function for the MC sampling, and includes the phase information in the accumulation. 

In Ref. H a general scheme for tackling the sign problem in PIMC simulations was proposed. The method has been 
applied to interacting electrons in a quantum dot Q and to the real-time dynamics of simple few-degrees-of-freedom 
- ■ systems Q . The generalization of the algorithm to the case of effective actions — which are long-ranged along the 
Trotter or real-time direction — along with the application to the dissipative two-state (spin-boson) dynamics has 
been given in Refs. (|]|]. This multilevel blocking (MLB) approach represents the systematic implementation of a 
simple blocking strategy. The blocking strategy states that by sampling groups of paths ("blocks") at the same time, 
the sign problem can always be reduced compared to sampling single paths as would be done normally; for a proof, 



see Sec. II A below. By suitably bunching paths together into sufficiently small blocks, the sign cancellations among 
paths within the same block can be accounted for without the sign problem, simply because there is no sign problem 
for a sufficiently small system. The MLB approach is based on a recursive implementation of this idea, i.e. after 
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building "elementary" blocks, new blocks are formed out of these, and the process is then iterated. This leads to a 
powerful implementation of the blocking strategy. This algorithm is able to turn the exponential severity of the sign 
problem into an "only" algebraic one. This is still difficult enough but in practice implies that significantly larger 
systems (lower temperature, longer real time) can be studied by MLB-PIMC than under the naive PIMC approach. 
Nevertheless, it should be stressed that MLB is definitely not a black-box scheme. There are several potential pitfalls 
related to incorrect or inadequate choices of certain MLB parameters, and one needs to be quite careful in applying 
this technique |10|. Given some experience, however, it represents a powerful handle to relieve the sign problem, with 
potential for additional improvement. 

The plan for the remainder of these notes is as follows. In Sec. || the MLB approach is discussed in detail, first 
in an intuitive way in Sec. [I A, and subsequently on a more formal or technical level in Sec. E B| for fermions. 
The modifications for real-time simulations are summarized in Sec. |II C This is then followed by a discussion of 
the performance and the accuracy of this algorithm in Sec. II D. In Sec. HE, the generalization to effective-action 



problems is outlined. To demonstrate the power of this approach, MLB results are presented in Sec. III. As an example 
for the fermionic sign prob lem, low-temperature simulation results for strongly correlated electrons in a 2D quantum 
dot are shown in Sec. [II A . The remainder of that section is the n concerned with the dynamical sign problem. After 
presenting result s for s everal simple model systems in Sec. [II B , the dynamics of the dissipative two-state system is 
discussed in Sec. Ill C. Finally, some conclusions can be found in Sec. IV. 



II. MULTILEVEL BLOCKING (MLB) APPROACH 



Before diving into the details of the MLB approach, the underlying idea ("blocking strategy") will be explained, 
focusing for simplicity on fermionic imaginary-time simulations. For those interested in working with this method, 
technical details and practical guidelines are provided in Sees. [IB to [ID. In the last part the generalization of MLB 
to PIMC simulations of the effective-action type is described. 



A. Blocking strategy 

Let us consider a many-fermion system whose state is described by a set of quantum numbers f denoting, e.g. the 
positions and spins of all particles. These quantum numbers can correspond to electrons living on a lattice or in 
continuous space. For notational simplicity, we focus on calculating the equilibrium expectation value of a diagonal 
operator or correlation function (this can be easily generalized), 

where Y^-. represents either a summation for the case of a discrete system or an integration for a continuous system, 
and p(r, r*) denotes the (reduced) density matrix of the system. In PIMC applications, imaginary time is discretized 
into P slices of length e = (3/ P. Inserting complete sets at each slice m = 1, . . . , P, and denoting the corresponding 
configuration on slice m by f m , the diagonal elements of the density matrix at f — fp entering Eq. (|l|) are: 

p 

p(rp,r P )= Y; l[(r m +i\e- eH \r m ) . (2) 
f t x,...,rp—i m=X 

Of course, periodic boundary conditions have to be enforced here, fp + i = r\. To proceed, one then has to construct 
accurate analytical approximations for the short-time propagator. This formulation of the problem excludes effective 
actions such as those arising from an integration over the fermions using the Hubbard-Stratonovich transformation 
H, since that generally leads to long-ranged imaginary-time interactions. The MLB approach suitable for such a 



situation l8f| is described below in Sec. [I E 



Since we are dealing with a many-fermion system, the short-time propagators need to be antisymmetrized, leading 
to the appearance of determinants causing the sign problem. Strictly speaking, the antisymmetrization has to be 
done only on one time slice, but the intrinsic sign problem is much better behaved if one antisymmetrizes on all time 
slices. Choosing the absolute value of the product of the short-time propagators in Eq. (|^) as the positive definite MC 
weight function one has to keep the sign Q[X] associated with a particular discretized path X = (r*i, . . . ,rp) 

during the accumulation procedure, 
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(A) 



J2 X p[x]$[x] 



(3) 



Assuming that there are no further exclusivity problems in the numerator so that A[X] is well-behaved, one can then 
analyze the sign problem in terms of the variance of the denominator, 



1 

~N S 



«$ 2 



(4) 



where N s is the number of MC samples taken and stochastic averages are calculated with P as the weight function. 
For the fermion sign problem, where $ = ±1 and hence ($ 2 ) = 1, the variance of the signal is controlled by the size 

of K*>|. 

Remarkably, one can achieve considerable progress by blocking paths together. Instead of sampling single paths 
along the MC trajectory one can consider sampling sets of paths called blocks. Under such a blocking operation, the 
stochastic estimate for (A) takes the form 



(A) 



J2b (Lxe B P[XMX]A[X]) 
\J2xeB 



(5) 



where one first sums over the configurations belonging to a block B in a way that is not affected by the sign problem, 
and then stochastically sums over the blocks. The summation within a block must therefore be done non-stochastically, 
or alternatively the block size must be chosen sufficiently small. Of course, there is considerable freedom in how to 
choose this blocking. 

Let us analyze the variance a' 2 of the denominator of Eq. (^|). First define new sampling functions in terms of the 
blocks which are then sampled stochastically 



P'[B] 



E p\xmx] 



xeB 



$'[B]= S gn ^P[I]$[I] 

\X£B J 



(6) 



Rewriting the average sign in the new representation, i.e. using P'[B] as the weight, then inserting the definition of 
P' and $' in the numerator, 

wim = ^bP'\bW\b] = E X P[XMX] 
[ 1 u J2bP'[b] EbP'[b] ' 

and comparing to the average sign in the standard representation using P[X] as the weight, one finds 

l(<f')l _ J2 x p[x] 



m\ z b p'[b] 



By virtue of the Schwarz inequality 



E^ = E 



E P{XMX] 



xeB 



< 



E 



e™ 



xeB 



= E p w 



X 



it follows that for any kind of blocking, the average sign improves, > \(&)\- Furthermore, since (<£>' 2 ) = (<£> 2 } = 1, 

Eq. (0) implies that 



a < a 



(7) 



and hence the signal-to-noise ratio is always improved upon blocking configurations together. Clearly, the worst 
blocking one could possibly choose would be to group the configurations into two separate blocks, one collecting all 
paths with positive sign and the other with negative sign. In this case, blocking yields no improvement whatsoever, 
and the "<" becomes "=" in Eq. (0). It is apparent from Eq. ^ that the blocking strategy provides a systematic 
handle to reduce the sign problem. In our realization of the blocking strategy, a block is defined by all paths that 
differ only at one time slice, i.e. only r m is updated with all other r n ^ m being frozen. 

A direct implementation of the blocking strategy does indeed improve the sign problem but will not remove its 
exponential severity. The reason is simply that for a sufficiently large system, there will be too many blocks, and 
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once the signals coming from these blocks are allowed to interfere, one again runs into the sign problem, albeit with a 
smaller scale Eq entering the signal-to-noise ratio. The resolution to this problem comes from the multilevel blocking 
(MLB) approach |3| where one applies the blocking strategy in a recursive manner to the blocks again. In a sense, 
new blocks containing a sufficiently small number of elementary ones are formed, and this process is repeated until 
only one block is left. Each step of this hierarchy is called level in what follows. Blocks are then defined living on 
these levels, and after taking care of the sign cancellations within all blocks on a given (fine) level, the resulting sign 
information is transferred to the next (coarser) level. On each step, the blocking strategy ensures that no sign problem 
occurs provided one has chosen sufficiently small block sizes. By doing this recursively, the sign problem on all the 
coarser levels can be handled in the same manner. It is then possible to proceed without numerical instabilities from 
the bottom up to the top level. The cancellations arising at the top level will create a sign problem again, which 
is however strongly reduced. As is argued below, the resulting sign problem is characterized by an only algebraic 
severity. 

In many ways, the MLB idea is related to the renormalization group approach. But instead of integrating out 
information on fine levels, sign cancellations are "synthesized" within a given level and subsequently their effects are 
transferred to coarser levels. While the renormalization group filters out information judged "relevant" and then 
ignores the "irrelevant" part, no such approximation is introduced here. Therefore our approach is actually closer in 
spirit to the multi-grid algorithm jll[| . The technical implementation of MLB is discussed next following Ref. ||] . 




FIG. 1. Levels for L = 2 (P = 4). Imaginary time flows along the circle (solid curve), and the slices m = 1,2,3,4 are 
distributed among the three levels: The finest level £ — contains m = 1,3, level 1=1 contains m = 2, and £ = 2 contains 
m = 4. Levels bonds are indicated by dashed and dotted lines. 

B. Systematic implementation: MLB 

To keep notation simple, the slice index m is used as a shorthand notation for the quantum numbers r m . From 
Eq. (^) the level-0 bonds, which are simply the short-time propagators, then follow in the form: 

(m,m + l) = (r m+1 \e- eH \r m ) . (8) 

Next the different levels < I < L, where L defines the Trotter number P = 2 L , have to be specified. Each 
slice m belongs to a unique level I, such that m = (2j + 1)2 £ and j is a nonnegative integer. For instance, the 
slices m = 1, 3, 5, • • • , P — 1 belong to I = 0, m = 2, 6, 10, • • • , P — 2 belong to I = 1, etc., such that there are 
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Mt = 2 L ~ l ~ 1 (but Ml — 1) different slices on level £, see Fig. |[ An elementary blocking is achieved by grouping 
together configurations that differ only at slice m, so only r m varies in that block while all r m /^ m remain fixed. 
Sampling on level £ therefore extends over configurations {r m } living on the Me different slices. In the MLB scheme, 
one moves recursively from the finest (£ = 0) up to the coarsest level (£ = L), and the measurement of the diagonal 
operator is done only at the top level using the configuration rp. 

A MC sweep starts by changing only configurations associated with the slices on level I — according to the 
standard weight 

Pb = |(l,2) (2,3)o---(P,l)o|, (9) 

generating a MC trajectory containing K samples for each slice on level £ = 0. These MqK samples are stored and 
they are used to generate additional coarser interactions among the higher-level slices, 

(m, m + 2)i = (sgn[(m, m + l) (m + 1, to + 2) Q ]) Vo[m+1] = K^ 1 ^ sgn[(m, to + l)o(m + 1, m + 2) ] , (10) 

m+l 

where the summation J2 m +i extends over the K samples available for slice to. As will be discussed later on, the 
important MLB parameter K — subsequently called sample number — should be chosen as large as possible to ensure 
that the second equality in Eq. ( |To| ) is justified. The level- 1 bonds ( |Io| ) contain precious and crucial information about 
the sign cancellations that occured on the previous level t = 0. Using these bonds, the density matrix (||) is rewritten 
as 

P (P,P)= J2 l(M)o(2,3)o ■•• (P,l)o\(2A)i •••(P-2, P)i(P,2)i . (11) 

1,2,...,F-1 

Comparing this to Eq. (||), the entire sign problem has been transferred to the next coarser level by using the level-1 
bonds. 

In the next step, the sampling is carried out on level £ — 1 in order to generate the next-level bonds, i.e. only slices 
m = 2, 6, . . . , P — 2 are updated, using the weight VoPi with 

Vi = |(2,4) 1 (4,6)i---(P,2) 1 | . (12) 

Moving the level-1 configurations modifies the level-0 bonds, which in turn requires that the level-1 bonds be updated. 
A direct re-calculation of these bonds according to Eq. (|l(]) would be too costly. Instead, the stored configurations 
on level £ = are used to perform an importance sampling of the new level-1 bonds. Under the test move m — > m', 
i.e. r m — > r' m , on level £ = 1, the bond ( |To| ) can be obtained from 

(m / ,m+l) (m+l,m+2)n 

<™' rr, J- 9\ — ^ m+1 l(m,m+l)o(m+l,m+2) | „s 
{m,m+*)l-— |( m ', m+ l) (m+l, m+ 2) | ' ^) 

£-im+\ | (m,m+l)o(m+l,m+2) 1 

where J2 m +i runs over the previously stored MC configurations r m+1 . Note that for small values of K, Eq. ( |l3| ) is 
only approximative, and thus a sufficiently large value of K should be chosen. With the aid of Eq. (13), the updated 
level-1 bonds follow with only moderate computational effort. Generating a sequence of K samples for each slice on 
level £ = 1 and storing them, level-2 bonds can be calculated in analogy to Eq. (|l0|), 

(m,m + 4) 2 = (sgn [(to, to + 2)i(m + 2, to + A)i]) VlVo . (14) 

Finally, the process is iterated up to the top level £ — L using the obvious recursive generalization of Eqs. ([l0|) and 
( |l4| ) to define level-£ bonds. 

Thereby the diagonal elements of the density matrix are obtained as 

P {P,P)= l(l,2)o(2,3)o---(F,l) | |(2,4)i...(P-2,P) 1 (P,2) 1 |...|(P/2,P) L _ 1 (P,P/2) i _ 1 |(P,P) i . (15) 

1,2,...,P-1 

By virtue of this algorithm, the sign problem is transferred step by step up to the coarsest level. The expectation 
value (0) can thus be computed from 

{A(P)sgn(P,P) L )v 
{A) (sgn(P,P) L ) P ' (16) 
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The manifestly positive definite MC weight V used for the averaging in Eq. ([l6]) can be read off from Eq. ( fl5| ) , 

V = |(1, 2) (2, 3) •• • (P, l)o] 1(2, 4)i • • • (P - 2, P)i(P, 2)x| • ■ ■ |(P/2, P)l-i(P, P/2)l-i\ |(P, P)l| • (17) 

The denominator in Eq. (|l6| ) gives the average sign and indicates to what extent the sign problem has been solved. 
For proper choice of MLB parameters, in particular the sample number this method can solve the sign problem. 
The price to pay for the stability of the algorithm is the increased memory requirement ~ K 2 associated with having 
to store the sampled configurations. 

C. Real-time simulations 

The same method, described so far for fermions, can be applied with minor modifications to a computation of real- 
time correlation functions or occupation probabilities. For concrctcness, let us focus on an equilibrium time-correlation 
function, 

( A (°)B(t)) = ^ {e _ lPh+a)H/he+aH/A} J ■ (18) 
Similarly other dynamical properties like the thermally symmetrized correlation function Jl2[|, 

c s (t) = z-^e^m^tw/h^-ip/m^tw/h^ ^ (19) 

with Z being the partition function, can be computed. In terms of path integrals, the traces in ( |l8| ) involve two 
quantum paths, one propagated backward in time for the duration — t and the other propagated in complex time for 
the duration t — i(3h. Discretizing each of the two paths into P slices, the entire cyclic path has a total of 2P slices. 
A slice on the first half of them has length —t/P and on the second half (t — ij3Ti)/ P. Denoting the quantum numbers 
(e.g. spin or position variables) at slice j by Tj, the correlation function ( |l8| ) reads 

J dri ■ ■ ■ dr 2P B(r 2 p)A{rp) ]jg ifa, r j+1 ) 
fdrf- dr 2P JJ 2 ^ (r, , r j+1 )„ 

where the level-0 bond (rj,r J+ i)o is again the short-time propagator between slices j and j + 1, and r 2 p+\ = r\. 
First assign all slices along the discretized path to different levels I = 0, . . . , L, where P — 2 L , in close analogy to the 
treatment for fermions, see Fig. [j]. Each slice j — 1, . . . , 2P belongs to a unique level £, such that j = (2k + l)2 e and 
k is a nonnegative integer. For instance, slices j = 1, 3, 5, . . . belong to level £ — 0, slices j = 2, 6, 10, . . . to I = 1, 
etc. The MLB algorithm starts by sampling only configurations which are allowed to vary on slices associated with 
the finest level £ = 0, using the weight Vq = |(ri, 7*2)0 ■ ■ ■ ( r 2P, r i)o|- The short-time level-0 bonds are then employed 
to synthesize longer-time level-1 bonds that connect the even-j slices. Subsequently the level-1 bonds are used to 
synthesize level-2 bonds, and so on. In this way the MLB algorithm moves recursively from the finest level (£ = 0) 
up to increasingly coarser levels until £ = L, where the measurement is done using r 2 p and rp. 

Generating a MC trajectory containing K samples for each slice on level I = and storing these samples, the level-1 
bonds © can be computed, where the "sgn" has to be replaced by the complex- valued phase <i>[z] = e mrg ^'. Their 
benefit becomes clear when rewriting the integrand of the denominator in ( po|) as 

Po x [r 2 ,r 4 ) 1 ■ ■ ■ (r 2P - 2 ,r 2 p) 1 (r 2P ,r 2 ) 1 . 

Comparing this to (pp[), the entire sign problem has been transferred to the next coarser level. In the next step, 
the sampling is carried out on level £ — 1 in order to compute the next-level bonds, using the weight P0P1 with 
Pi = |( r 2,'"4)i • • • {r 2 p, r 2 )i \. Generating a sequence of K samples for each slice on level 1=1, and storing these 
samples, level-2 bonds can be calculated, 

{rj,r j+i ) 2 = ($[(^,^+2)1(^+2, Tj+^iDvoV! ■ 
This process is then iterated up to the top level. Finally, the correlation function ( |T8| ) can be computed from 

(B(r 2P )A(rp)${(rp,r 2 p) L (r 2 p,rp) L ])-p 



($[(r P , r 2P ) L (r 2P , r P ) L ]) v 



(21) 



with the positive definite MC weight P = P0P1 • ■ -Pl- The denominator in Eq. ( pl| ) gives the average phase and 
indicates to what extent the sign problem has been solved. With a suitable choice of MLB parameters, the average 
phase remains close to unity even for long times. 
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D. Accuracy and pitfalls 



Next questions concerning the exactness and performance of the MLB approach are addressed, in particular the 
dependence on the sample number K. Clearly, K needs to be sufficiently large to produce a reliable estimate for the 
level- ^ bonds. If these bonds could be calculated exactly — corresponding to the limit K — > oo — , the manipulations 
leading to Eq. yield the exact result. Hence for large enough K, the MLB algorithm must become exact and 
completely solve the sign problem. Since the levels bonds can however only be computed for finite K, the weight 
function V amounts to using a noisy estimator, which in turn can introduce bias into the algorithm p3[. In principle, 
this problem could be avoided by using a linear acceptance criterion instead of the algorithmically simpler Metropolis 
choice Q which was employed in the applications reported here. But even with the Metropolis choice, the bias can be 
made arbitrarily small by increasing K . Therefore, with sufficient computer memory, the MLB approach can be made 
to give numerically exact results. One might then worry about the actual value of K > K* required to obtain stable 
and exact results. If the critical value K* were to scale exponentially with (3 and/or system size, the sign problem 
would be present in disguise again. 

Although a rigorous non-exponential bound on K* has not yet been established, our experience with the MLB 
algorithm indicates that this scaling is at worst algebraic. This is corroborated by a recent careful study of this issue 
|[10|| . Therefore the exponential severity of the sign problem is replaced by an algebraic one under MLB. A heuristic 
argument supporting this statement goes as follows. If one needs K samples for each slice on a given level in order 
to have satisfactory statistics despite of the sign problem, the total number of paths needed in the naive approach 
depends exponentially on P, namely ~ K p . This is precisely the well-known exponential severity of the sign problem 
under the naive approach. However, with MLB the work on the last level, which is the only one affected by a sign 
problem provided K was chosen sufficiently large, scales only ~ K L . Note that it does not scale ~ K because one 
must update the level- £ bonds on all L finer levels as well. So in MLB, the work needed to sample the K p paths with 
satisfactory statistical accuracy grows ~ K lo&2 p = P lo S2 K t j. e . algebraically with P. An important point to mention 
at this point concerns the high-temperature (or short-time) limit, where the direct PIMC simulation does not face 
a significant sign problem. In this case, however, the above-mentioned bias problematics of the MLB-PIMC is quite 
serious and can give erroneous results. Fortunately, since that regime is of little interest to MLB, this is not a serious 
restriction. A more detailed discussion of this point can be found in Ref. (h}. 

To elucidate how the MLB algorithm works in practice, in Table |l| sim ulation results for N — 8 electrons in a 



quantum dot at various values of K are listed. For details, see Sec. [II A. Compared to the naive approach where 
K = 1, using a moderate K — 200 already increases the average sign from 0.02 to 0.63, making it possible to obtain 
more accurate results from much fewer samples. The data in Table | also confirms that the bias can be systematically 
eliminated by increasing K , so that the energy found at K > 200 essentially represents the exact result (within error 
bars). The value K* « 200 is quite typical for many applications. For a simple model system, a value K* ps 50 
was found in Ref. |10|. Table [j] also shows that the CPU time per sample scales linearly with K, whereas memory 
requirements grow ~ K 2 . 



TABLE I. MLB results for N = 8 and A = 2 (see Sec. |lIAj). N s is the number of samples (in 10 4 ), t C p\j the total CPU 
time (in hours), MB the required memory (in mega-bytes), and (sgn) the average sign. Bracketed numbers are error estimates. 



K 
1 

100 
200 
400 
600 
800 



N 3 

120 

7 

9 

8 

10 



£cpu 

95 

33 

83 

174 

308 

429 



MB 

2 

14 
30 
64 
96 
150 



(sgn) 

0.02 

0.48 

0.63 

0.73 

0.77 

0.81 



En /hup 

48.6(3) 

48.43(8) 

48.55(7) 

48.53(9) 

48.54(8) 

48.59(8) 
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E. Effective actions 



So far the MLB algorithm was only discussed for the case of nearest- neighbor interactions along the Trotter/time 
direction. This situation is encountered under a direct Trotter-Suzuki breakup of the short-time propagator. In many 
cases, however, one has to deal with effective actions that may include long-ranged interactions along the (complex) 
time direction. Such effective actions arise from degrees of freedoms having been traced out, e.g. a harmonic heat 
bath |l4|], or through a Hubbard-Stratonovich transformation in auxiliary-field MC simulations of lattice fermions 
PJ. Remarkably, because such effective actions capture much of the physics such as symmetries or the dissipative 
influence of traced-out degrees of freedom, the corresponding path integral very often exhibits a significantly reduced 
intrinsic sign problem compared to the original (time-local) formulation. To be specific, let us focus on the dynamical 
sign problem arising in real-time PIMC computations here. The modifications required to implement the method for 
fermion simulations are then straightforward. 

Let us consider a disc retize d path integral along a slightly different but fully equivalent contour in the complex-time 



plane compared to Sec. [I C , namely a forward branch from t = to t = t*, where t* is the maximum time studied 
in the simulation, followed by a branch going back to the origin, and then by an imaginary-time branch from t = 
to t = —i%j3. Here a "factorized" initial preparation is studied, where the relevant degrees of freedom, r(i), are held 
fixed for t < Jf^ ]. That implies that the imaginary-time dynamics must be frozen at the corresponding value, and 
one only needs to sample on the two real-time branches. Note that such a nonequilibrium calculation cannot proceed 
by first doing an imaginary-time QMC simulation followed by a generally troublesome analytic continuation of the 
numerical data 0. Using time slices of length t*/P, forward [r(i m )] and backward [r'(t m )) path configurations at 
time t rn — mt* / P are combined into the configuration s m , where m = 1, . . . , P. The configuration at t = is held 
fixed, and for t = t* one must be in a diagonal state, r(t*) = r'(t*). For an efficient application of the current method, 
it is essential to combine several neighboring slices m into new blocks. For instance, think of m — 1, . . . , 5 as a new 
"slice" I = 1, m = 6, . . . , 10 as another slice I = 2, and so on. Combining q elementary slices into a block sg, instead 
of the original P slices one has L = P/q blocks, where L is the number of MLB levels. Instead of the "circular" 
structure of the time contour inherent in the trace operation, it is actually more helpful to view the problem as a 
linear chain, where the MLB scheme proceeds from left to right. In actual applications, there is considerable freedom 
in how these blocks are defined, e.g. if there is hardly any intrinsic sign problem, or if there are only few variables in 
r, one may choose larger values of q. Additional flexibility can be gained by choosing different q for different blocks. 

Say one is interested in sampling the configurations sl on the top level £ = L according to the appropriate matrix 
elements of the (reduced) density matrix, 

p(s L ) = Z- 1 exp{-5[fli,...,s L ]} , (22) 

Sl,...,3L-l 

where S is the effective action under study and Z is a normalization constant so that ^2 S p(sl) = 1- Due to the 
time- non-locality of this action, there will be interactions among all blocks sg. The sum in Eq. ( p2| ) denotes either 
an integration over continuous degrees of freedom or a discrete sum. In the case of interest here, the effective action 
is complex- valued and e~ s /\e~ s \ represents an oscillatory phase factor $. The "naive approach" to the sign problem 
is to sample configurations using the positive definite weight function V ~ | exp{— S}\, and to include the oscillatory 
phase in the accumulation procedure. Below it is assumed that one can decompose the effective action according to 

L 

S[s 1 ,...,s L ]=Y,W i [se,...,s L ] . (23) 
t=i 

All dependence on a configuration si is then contained in the "partial actions" W\ with A < £. One could, of course, 
put all Wf>i = 0, but the approach becomes more powerful if a nontrivial decomposition is possible. 

Let us now describe the algorithm in some detail, employing a somewhat different but equivalent notation than 
before. This may be helpful to some readers in order to better understand the MLB algorithm, see also Ref. 0] 
for a formulation of Sec. [II B] in this notation. The MC sampling starts on the finest level I = 1, where only the 
configuration si = i containing the elementary slices m — 1, . . . , q will be updated with all s^>i remaining fixed at their 
initial values s°. Using the weight function 



V [si] = \exp{-Wx[s 1 ,s 2 ,...,s° L \}\ 

. . . , K, and store them for later use. 
be chosen large enough. For K = 1, the algorithm will simply reproduce the naive approach. The stored samples 



generate K samples , where i = 1, . . . ,K, and store them for later use. As usual, the sample number K should 
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are now employed to generate information about the sign cancellations. All knowledge about the interference that 
occured at this level is encapsulated in the quantity 



B, = 



cxp{-Wi[si, . . . , s L }} 
exp{-W 1 [s 1 ,s 2 ,...,s L ]}\ 



Vo[si 



= c o ^^^expi-W^Si, . . . , s L }} 



(24) 



A" 



exp{-Wi[s^ , s 2 , . . . , s L ]} 
~[ \ exp{-W 1 [s ( i\s 2 ,..., S l}}\ 



= B x [s 2 , ...,s L ] , 



which are analogously called level-1 bonds, with the normalization constant Co = ^ot^i]- Combining the second 
expression in Eq. (E4T) with Eq. (p3), the density matrix reads 



P (s L ) = Z- 1 ]T expl-J^wAcaBx = Z' 1 £ V^R, 

S2,-,sl-i I £>1 ) si,...,ai_i 1>1 



(25) 



When comparing Eq. ( |25| ) with Eq. fl22|), the sign problem has now been transferred to levels I > 1, since oscillatory 
phase factors only arise when sampling on these higher levels. Note that B\ = Bi[s 2 , ■ ■ ■ , s_l] introduces couplings 
among all levels I > 1, in addition to the ones already contained in the effective action S. 

Next proceed to the next level I — 2 and, according to Eq. (p5|), update configurations for m = q + 1, . . . , 2q using 
the weight 



(26) 



Under the move s 2 — > s' 2 , one should then resample and update the level-1 bonds, B\ — » B[. Exploiting the fact 
that the stored K samples sf* are correctly distributed for the original configuration s 2 , the updated bond can be 
computed according to 



A" 



B[=K-^ 



exp{-W 1 [s [ *\s' 2 ,...,s L )} 
t \eM-W l [sf\s° 2 ,..., S l}}\ 



(27) 



Again, to obtain an accurate estimate for BJ, the number K should be sufficiently large. In the end, sampling under 
the weight V\ implies that the probability for accepting the move s 2 —* s' 2 under the Metropolis algorithm is 



V 



cxp{-W 1 [ 8 W,^, S °,...]} 



cxp{-Wi[si 



]}l 



cxp{-Wi[ar,82,ag,-]} 



cxp{~W 2 [s' 2 ,sl...}} 



cxp{-W 2 [s 2 ,s° 3 ,...]} 



(28) 



'* |exp{-W 1 [*i i, , a S>,...]}| 

Using this method, one generates samples 3 stores them, and computes the level-2 bonds, 

j-h = i r;:r 1 : ::r~-~:;'" ; > =cr i ^s 1 [ S2) ...]ex P {-w 2 [ S2) ...]} 



gijggvfs, ■ ■ ■] exp{-VF 2 [s 2 , g 3, ■ ■ ■]} \ 
iBifaa, •§,...] exp{-^ 2 [s 2 , S o,...]}|/ Pi[s2 



(29) 



" j?i[4 t) ^3,...]exp{-lU 2 [4 i) ,8 3 ,...]} 
|Sa[«W, «g, . - .] exp{-Wa[^, «g, . . .]>j 



= B 2 [s 3) ...,s L ] 



with C\ = ^2 S2 "Pi[s 2 ]- Following above strategy, one then rewrites the reduced density matrix by combining Eq. fl25| ) 
and the second expression in Eq. (129), 



p(s L ) = z-i J2 



exp ■ 



S3,...,SL-1 



£>2 



W t > C CiB 2 = Z- 



E 

ei,...,BL—i 



w, 



(30) 



e>2 



Clearly, the sign problem has been transferred one block further to the right along the chain. Note that the nor- 
malization constants Co, Cl, . . . depend only on the initial configuration s° so that their precise values need not be 
known. This procedure is then iterated in a recursive manner. Sampling on level I using the weight function 



Ve-i[se] = \Bt-i[st, s° e+1 , . . .} exp{-W e [s e , s° +1 , . . .]} 



(31) 
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requires the recursive update of all bonds B\ with A < £. Starting with B\ — > B[ and putting Bq — 1, this recursive 
update is done according to 



B , _ K -i y- B' x _ 1 [8 { £ ) ,a x+1 , . . .} exp{-W{[s ( £ ) ,a x+1 , . . .}} 
>X " U \B x . 1 [s ( ^\sl +1: ...}eM-W x [s^,sl +17 ...}}\ 1 



where the primed bonds or partial actions depend on s' t and the unprimed ones on a®. Iterating this to get the 
updated bonds Bi_ 2 for all a^_ 1} the test move si — > s' t is then accepted or rejected according to the probability 



B t -i[a' t , 8° e+1 , . . .] exp{-H^[4, s° +1 , . . .]} 



B t -i[s t , s° e+1 , . . .] exp{-Wt[at, ■ ■ ■]} 
On this level, one again generates K samples sf\ stores them and computes the level-£ bonds according to 



(33) 



Bi[s e+1 , ...}= K 



_i Bi^i[s\ L> ,si + i, . . .}exp{-Wi[s ( t\s e+1 , . . .]} 
U \B t ^\ 8 ° t+1 , . . .} exp{-Wrf, ■ ■ .]}| 



This process is iterated up to the top level, where the observables of interest may be computed. Since the sampling 
of Bi requires the resampling of all lower-level bonds, the memory and CPU requirements of the algorithm laid out 
here are quite large. For A < £ — 1, one needs to update B\ — > B' x for all si, with A < £' < t, which implies a 
tremendous amount of computer memory and CPU time, scaling approximately ~ K L at the top level. Fortunately, 
an enormous simplification can often be achieved by exploiting the fact that the interactions among distant slices 
are usually weaker than between near-by slices. For instance, when updating level t — 3, the correlations with the 

(i) . . (i) 

configurations s\ may be very weak, and instead of summing over all K samples s\ in the update of the bonds 
Ba<£, one may select only a small subset. When invoking this argument, one should be careful to also check that the 
additional interactions coming from the lcvcl-A bonds with A < I are sufficiently short-ranged. From the definition of 
these bonds, this is to be expected though. 



III. APPLICATIONS 



In this section, several different applications of the MLB approach will be presented. The first will focus on the 
equilibrium behavior of interacting electrons in a parabolic quantum dot, a situation characterized by a fermionic sign 
problem. The two other subsections then deal with the dynamical sign problem. 



A. Quantum dots 

Quantum dots are solid-state artificial atoms with tunable properties. Confining a small number of electrons N 
in a 2D electron gas in semiconductor heterostructures, novel effects due to the interplay between confinement and 
the Coulomb interaction have been observed experimentally []ili|-|l7f . For small N, comparison of experiments to the 
generalized Kohn theorem indicates that the confinement potential is parabolic and hence quite shallow compared 
to conventional atoms. Employing the standard electron gas parameter r s to quantify the correlation strength, only 
for small r s , a Fermi- liquid picture is applicable. In the low-density (strong-interaction) limit of large r s , classical 
considerations suggest a Wigner crystal-like phase with electrons spatially arranged in shells. We call this a Wigner 
molecule due to its finite extent. Of particular interest is the crossover regime between these two limits, where both 
single-particle and classical descriptions break down, and basically no other sufficiently accurate method besides QMC 
is available. Exact diagonalization is limited to very small N since one otherwise introduces a huge error due to the 
truncation of the Hilbert space. Hartree-Fock (and related) calculations become unreliable for large r s and incorrectly 
favor spin-polarized states. Furthermore, density functional calculations can introduce uncontrolled approximations. 
Regarding QMC, to avoid the sign problem, the fixed-node approximation has been employed by Bolton Jig ] and 
later by others |l^]. For N > 5, typical fixed-node errors in the total energy are found to be of the order of 10% 
||. It is then clear that one should resort to exact methods, especially when looking at spin-dependent quantities, 
where often extremely small spin splittings are found. After our original studies |]J(|, other studies using the naive 
PIMC approach were published p0pl|] . These studies are however concerned with the deep Wigner regime r s > 20 
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pj| , which is essentially a classical regime without sign problem not further discussed here, or employ a special 
virial estimator [^0) that unfortunately appears to be incorrect except for fully spin-polarized states p2[ . A clean 2D 
parabolic quantum dot in zero magnetic field is described by 

The electron positions (momenta) are given by Xj (pj), their effective mass is to*, and the dielectric constant is k. 
The MLB calculations are carried out at fixed N and fixed z-component of the total spin, 5" = (A-f — AjJ/2. As a 
check, the exact solution for A = 2 has been reproduced. 

Here results for the energy, E — (H) (since H is a non-diagonal operator, two Trotter slices are kept at the top 
level), and the spin sensitivity ^{r s ) will be discussed. The latter quantity is useful to study the crossover from weak 
to strong correlations, 



/>oo 

i N (r s ) oc y / dyylas 

s,s> Jo 



(y)-g S '(y)\, (35) 



where the prefactor is chosen to give £at = 1 for r s — 0. This definition makes use of the spin-dependent two-particle 
correlation function 



9s( x ) 



2-Kll 



N V 

y S(x-xt + Xj) J , (36) 



N{N - 1) 

which is isotropic. With y = r/lo prefactors are chosen such that dyyg s (y) = 1. The confinement length scale Z = 
^/h/m*ojQ allows the interaction strength to be parametrized by the dimensionless parameter A = la /a = e 2 /K,u>olo, 
where a is the effective Bohr radius of the artificial atom. For any given N and A, the density parameter r s ~ r* /a 
with nearest-neighbor distance r* follows by identifying r* with the first maximum in 9s( r )- The correlation 
function ( |36"|) is a very sensitive measure of Fermi statistics, in particular revealing the spin-dependent correlation 
hole. Because interactions tend to wash out the Fermi surface, the spin sensitivity ( |35| ) is largest for a Fermi gas, 
r s = 0. Since for r s — > oo, one approaches the totally classical limit, where gs (r) is completely spin- independent, 
^n(ts) decays from unity at r s = down to zero as r s — > oo. The functional dependence of this decay provides insight 
about the crossover phenomenon under study. 
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FIG. 2. Numerical results for £jv(r s ) at hsT/huc — 0.1. Statistical errors are of the order of the symbol size. The dotted 
curve, given by exp(— r s /r c ) with r c — 4, is a guide to the eye only. The inset shows the same data on a semi-logarithmic scale 
as a function of ^/rj. The dashed line is the WKB estimate (see text). 

By computing charge densities, the PIMC simulations can directly reveal shell formation in real space ||. Such a 
spatial structure is clear evidence for near-classical Wigner molecule behavior. The classical shell filling sequence is 
as follows. For N < 6, the electrons arrange on a ring, but the sixth electron then goes into the center. Furthermore, 
electrons 7 and 8 enter the outer ring again. These predictions are in accordance with our PIMC data. Clear indications 
of a spatial shell structure at N > 6 can be observed already for A « 4, albeit quantum fluctuations tend to wash 
them out somewhat. For A>4, charge densities are basically insensitive to S. This is characteristic for a classical 
Wigner crystal, where the Pauli principle and spin-dependent properties are of little importance. Numerical results 
for the spin density in this regime simply follow the corresponding charge density according to s z (r) ~ (S/N)p(r). 
A significant S'-dependence of charge and spin densities is observed only for weak correlations. Figure || reveals that 
£jv( r s) i s remarkably universal, i.e. it depends only very weakly on N. Its decay defines a crossover scale r c , where 
an exponential fit for small r s yields r c ks 4. For r s > 4, the data can be well fitted by £(r s ) ~ exp(— \/r s /r' c ), 
where r' c w 1.2. Remarkably, this is precisely the behavior expected from a semiclassical WKB estimate for a Wigner 
molecule |J . The crossover value r c « 4 is also consistent with the onset of spatial shell structures in the density, and 
with the spin-dependent ground state energies expected for a Wigner molecule. Therefore the crossover from weak to 
strong correlations is characterized by the surprisingly small value r c w 4, instead of r c w 37 found for the bulk 2D 
electron gas [^3| . This enormous stabilization of the Wigner molecule can be ascribed to the effects of the confinement 
potential. In the thermodynamic limit, loo — > with r s fixed, plasmons govern the low-energy physics, and hence the 
bulk value r c « 37 becomes relevant for very large N. For GaAs based quantum dots, one can estimate H that for 
N <, 10 4 , the value r c « 4 is the relevant one. Remarkably, recent experiments on vertical quantum dots jl(| have 
found evidence for an even smaller crossover scale r c w 1.8. The experimental study was carried out in a magnetic 
field, and the dot contained several impurities. Since both effects tend to stabilize a Wigner crystallized phase, our 
prediction and the experimental observation appear to be consistent. 
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TABLE II. MLB data for the energy for various {TV, S, A} parameter sets at k B T/hio = 0.1. Bracketed numbers denote 
statistical errors. 
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MLB results for the energy at different parameter sets {N, S, A} are listed in Table |J. For given N and A, if the 
ground state is (partially) spin-polarized with spin S, the simulations should yield the same energies for all S' < S. 
Within the accuracy of the calculation, this consistency check is indeed fulfilled. For strong correlations, r s > r c , the 
spin-dependent energy levels differ substantially from a single-particle orbital picture. In particular, the ground-state 
spin S can change and the relative energy of higher-spin states becomes much smaller than Tuoq. For N — 3 electrons, 
as r s is increased, a transition occurs from S = 1/2 to S = 3/2 at an interaction strength A ks 5 corresponding to 
r s ss 8. For N = 4, a Hund's rule case with a small-r s ground state characterized by S = 1 is encountered. From our 
data, this standard Hund's rule covers the full range of r s . A similar situation arises for N — 5, where the ground 
state is characterized by S = 1/2 for all r s . Turning to N — 6, while one has filled orbitals and hence a zero-spin 
ground state for weak correlations, for A = 8 a S = 1 ground state is found. A similar transition from aS = 1/2 state 
for weak correlations to a partially spin-polarized S = 5/2 state is found for N — 7. Finally, for N = 8, as expected 
from Hund's rule, a S = 1 ground state is observed for small r s . However, for A >4, corresponding to r s > 10, the 
ground state spin changes to S — 2, implying a different "strong-coupling" Hund's rule. 

Let us finally address the issue of magic numbers. For small r Sl the simple picture of filling up single-particle orbitals 
predicts that certain N are exceptionally stable. Results for the energy per electron, En/N, in the spin-polarized 
state S — N/2, are shown in Figure ||. Notably, there are no obvious cusps or breaks in the Af-dependence of the 
energy. The A = 2 data in Fig. || suggest that an explanation of the experimentally observed magic numbers [0 has 
to involve spin and/or magnetic field effects, since the single-particle picture breaks down so quickly. Remarkably, 
there are no pronounced cusps in E^/N for strong correlations (A = 8). Therefore magic numbers seem to play only 
a minor role in the Wigner molecule phase. 
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FIG. 3. Energy per electron, En/N, for S — N/2 and ksT/Tiujo — 1/6, in units of huo, for A = 2 (squares) and A = 8 
(diamonds). Statistical errors are smaller than the symbol size. Open circles are T — fixed-node QMC results Jl8| for A = 2. 
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B. Real-time simulations: Simple model systems 

In each of the following examples, a time-correlation function was computed directly in real time for a simple model 
system, with increasing level of complexity. The average phase is larger than 0.6 for the presented data sets. 

1. Harmonic oscillator 

For H = p 2 /2m + mw 2 ir 2 /2, the real and imaginary parts of (x(0)x(t)) oscillate in time due to vibrational coherence. 
In dimensionless units m — lo = 1, the oscillation period is 2tt. Figure |(a) shows MLB results for C(t) = Re (x(O)x(i)}. 
With P — 32 for the maximum time t = 26, K = 200 samples were used for sampling the coarser bonds. Within 
error bars, the data coincide with the exact result and the algorithm produces stable results free of the sign problem. 
Without MLB, the signal-to-noise ratio was practically zero for t > 2. 

2. Two-level system 

For a symmetric two-state system, H = —^Aa x , the dynamics is controlled by tunneling. The spin correlation 
function (a z (Q)a z (t)) exhibits oscillations indicative of quantum coherence. Figure f|(b) shows MLB results for C(t) = 
Re (<7z(0)<7z(i)), Putting A = 1, the tunneling oscillations have a period of 2tt. With P = 64 for the maximum time 
t = 64, only K = 100 samples were used for sampling the coarser bonds. The data agree well with the exact result. 
Again the simulation is stable and free of the sign problem. Without MBL, the simulation failed for t > 4. 

3. Double-well potential 

Next, let us examine a double-well system with the quartic potential V(x] = —x 2 + ^x 4 . At low temperatures, 
interwell transfer occurs through tunneling motions on top of intrawell vibrations. These two effects combine to 
produce nontrivial structures in the position correlation function. At high temperatures, interwell transfer can also 
occur by classical barrier crossings. Figure |J(c) shows MLB results for C(t) = Re (x(0)x(t)). The slow oscillation 
corresponds to interwell tunneling, with a period of approximately 16. The higher- frequency motions are characteristic 
of intrawell oscillations. In this simulation, K — 300 samples were used. The data reproduce the exact result well, 
capturing all the fine features of the oscillations. Again the calculation is stable and free of the sign problem, whereas 
a direct simulation failed for t > 3. 

4- Multidimensional tunneling system 

As a final example, consider a problem with three degrees of freedom, in which a particle in a double-well potential 
is bilinear ly coupled to two harmonic oscillators. The quartic potential in the last example is used for the double- well, 
and the harmonic potential in the first example is used for both oscillators. The coupling constant between each 
oscillator and the tunneling particle is a = 1/2 in dimensionless units. For this example, the correlation function 
C s (t) in Eq. jl9| ) has been computed for the position operator of the tunneling particle. Direct application of MC 
sampling to C s (t) has generally been found unstable for t > /3H/2 [jl2). In contrast, employing only moderate values 
of K = 400 to 900 allow to go up to t = 10/371. Figure |(d) shows MLB results for C' s {t) = ReC s {t). For the 
coupled system, the position correlations have lost the coherent oscillations and instead decay monotonically with 
time. Coupling to the medium clearly damps the coherence and tends to localize the tunneling particle. 
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FIG. 4. MLB results (closed circles) for various systems. Error bars indicate one standard deviation, (a) C(t) for a harmonic 
oscillator at [3Ti = 1. The exact result is indicated by the solid curve, (b) Same as (a) for a two-level system at f3Ti — 10. (c) 
Same as (a) for a double- well system at /3h = 1. This temperature corresponds to the classical barrier energy, (d) C' s (t) for a 
double- well system coupled to two oscillators at (5Ti — 1. For comparison, open diamonds are for the uncoupled (a = 0) system. 
Note that C' s (t) is similar but not identical to C(i) shown in (c). Solid and dashed lines are guides to the eye only. 



C. Spin-boson dynamics 

Finally, to demonstrate the performance of the MLB approach for effective-action-type problems, the real-time 
dynamics of the celebrated spin-boson model 



H = -(ftA/2) o x + (fie/2) a z + £ 



2m a + 2 a a 



(37) 



is discussed. This model has a number of important applications |14j], e.g. the Kondo problem, interstitial tunneling 
in solids, quantum computing and electron transfer reactions, to mention only a few. The bare two-level system 
(TLS) has a tunneling matrix element A and the asymmetry (bias) e between the two localized energy levels (a x 



1G 



and cr 2 are Pauli matrices). Dissipation is introduced via a linear heat bath, i.e. an arbitrary collection of harmonic 
oscillators {x a } bilinearly coupled to a z . Concerning the TLS dynamics, all information about the coupling to the 
bath is contained in the spectral density J{uj) — (n/2) J2 a ( c a/ m a (jJ a ) 5{uo — w a ), which has a quasi-continuous form 
in typical condensed-phase applications and dictates the form of the (twice-integrated) bath correlation function 

duj J(uj) cosh[w^/3/2] - cosh[cj(?i/?/2 - it)} 

^h~ 2 sinh[w?i/3/2] ' ( ' 

For the calculations here, an ohmic spectral density J(lu) = 2Trhau> exp(— uj/uj c ) has been assumed, for which Q(t) can 
be found in closed form [ fL4| . Here u> c is a cutoff frequency, and the damping strength is measured by the dimensionless 
Kondo parameter a. In the scaling limit A <C uj c with a < 1, all dependence on oj c enters via a renormalized tunnel 
splitting @ 

A eff = [cos(7ra)r(l - 2a)] 1/2(1 - Q) (AK)"^ 1 "^ A , (39) 

and powerful analytical and alternative numerical methods are readily available However, there are important 
applications (e.g. electron transfer reactions) that require to study the spin-boson problem away from the scaling 
limit. Here one generally has to resort to numerical methods. Basically all other available computational techniques 
can only deal with equilibrium quantities, or explicitly introduce approximations; for an overview and references, see 
Ref. [pR. The MLB approach is computationally more expensive than other methods but at the same time unique 
in yielding numerically exact results for the nonequilibrium spin-boson dynamics for arbitrary system parameters 
A,e,J(cj) and = l/k B T. 



TABLE III. 


MLB performance for a = 1/2, lj c /A = ( 


i, At* = 10, P = 40, 


and several L. 


qe denotes the number of slices 


for I = 1, . . . L. 










K 


L 


qe 




(sgn) 


1 


1 


40 




0.03 


200 


2 


30 


- 10 


0.14 


800 


2 


30 


- 10 


0.20 


200 


3 


22 


-12-6 


0.39 


600 


3 


22 


-12-6 


0.45 
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Below only results for the occupation probability P(t) = (o~ z (t)) under the nonequilibrium initial preparation o~ z (t < 
0) = +1 are presented. P(t) gives the time-dependent difference of the quantum-mechanical occupation probabilities 
of the left and right states, with the particle initially confined to the left state. To obtain P(t) numerically, in a 
discretized path-integral representation one traces out the bath to get a long-ranged effective action, the influence 



functional 14|. In discretized form the TLS path is represented by spins <7j, o~[ = ±1 on the forward- and backward- 
paths, respectively. The total action S consists of three terms. First, there is the "free" action So determined by the 
bare TLS propagator Uq, 



p-i 



exp(-So) = J] U (a l+1 ,af,t*/P) U Q (a' i+x , a'^ -t*/P) , (40) 

where t* is the maximum time and P the Trotter number. Next there is the influence functional, Si = Sj 1 ^ + Sj 2 \ 
which contains the long-ranged interaction among the spins, 

S\ 1] = ^2 ( a o ~ a 'j){ L 'j-m (°>n - o-' m ) + iL J _ m (a rn + a' m )j , 

j>m 

where Lj — L'j + iLj is given by 

Lj = [Q((j + l)t*/P) + Q((j - l)t*/P) - 2Q(ji*/P)]/4 (41) 

for j > 0, and Lo = Q(t* / P) / A. In the scaling regime at T = 0, this effective action produces interactions ~ a/t 2 
between the spins ("inverse-square Ising model"). The contribution 

sf = i(t*/P)J2l(mt*/P)(o-m -O 

m 

describes the interactions with the imaginary-time branch where a z = +1, with the damping kernel 

2 f°° , J(u) . . 
j{t) = — / duj cos(ujt) . 

The most difficult case for PIMC corresponds to an unbiased two-state system at zero temperature, e = T = 0. 
To check the code, the case a = 1/2 was studied in some detail, where the exact solution O] is very simple, 
P(t) = exp(— A c ffi). This exact solution only holds in the scaling limit, which is already reached for w c /A = 6 where 
the MLB-PIMC simulations yield precisely this result. Typical MLB parameters and the respective average sign are 



listed in Table III. The first line in Table [II corresponds to the naive approach. It is then clear that the average 
sign and hence the signal-to- noise ratio can be dramatically improved thus allowing for a study of long timescales t* . 
For a fixed number of levels L, the average sign grows by increasing the parameter K . Alternatively, for fixed K, the 
average sign increases with L. Evidently, the latter procedure is more efficient in curing the sign problem, but at the 
same time computationally expensive. In practice, it is then necessary to find a suitable compromise. 

Figure [5] shows scaling curves for P(t) at a = 1/4 for lj c / A = 6 and lj c /A = 1. The first value for w c /A is within 
the scaling regime. This is confirmed by a comparison to the noninteracting blip approximation (NIBA) JL4|, which 
is known to be very accurate for a < 1 in the scaling limit. However, for uj c /A = 1, scaling concepts and also NIBA 
are expected to fail dramatically. This is seen in the simulations. MLB results show that away from the scaling limit, 
quantum coherence is able to persist for much longer, and both frequency and decay rate of the oscillations differ 
significantly from the predictions of NIBA. In electron transfer reactions in the adiabatic-to-nonadiabatic crossover 
regime, such coherence effects can then strongly influence the low-temperature dynamics. One obvious and important 
consequence of these coherence effects is the breakdown of a rate description, implying that theories based on an 
imaginary-time formalism might not be appropriate in this regime. A detailed analysis of this crossover regime using 
MLB is currently in progress M. 
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A eff t 

FIG. 5. Scaling curves for P(t) for a — 1/4 with u> c /A — 6 (closed diamonds) and ui c /A = 1 (open circles). The solid curve 
is the NIBA prediction. Statistical errors are of the order of the symbol sizes. 

IV. CONCLUDING REMARKS 

These notes summarize our previous activities using the multilevel blocking approach to the sign problem in path- 
integral Monte Carlo simulations. The approach holds substantial promise towards relieving (and eventually over- 
coming) the sign problem, but clearly there is still much room for improvement. The applications presented here 
demonstrate unambiguously that there are general and powerful handles to relieve the sign problem, even though a 
problem characterized by an intrinsic sign problem is still much harder than one without. We hope that especially 
young researchers will be attracted to work on this subject themselves. 
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